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ABSTRACT 

When modeling astrophysical fluid flows, it is often appropriate to discard the 
canonical magnetohydrodynamic approximation thereby freeing the magnetic field 
to diffuse with respect to the bulk velocity field. As a consequence, however, the 
induction equation can become problematic to solve via standard explicit techniques. 
In particular, the Hall diffusion term admits fast-moving whistler waves which can 
impose a vanishing timestep limit. 

Within an explicit differencing framework, a multifluid scheme for weakly ionised 
plasmas is presented which relies upon a new approach to integrating the induction 
equation efficiently. The first component of this approach is a relatively unknown 
method of accelerating the integration of parabolic systems by enforcing stability over 
large compound timesteps rather than over each of the constituent substeps. This 
method, Super Time Stepping, proves to be very effective in applying a part of the 
Hall term up to a known critical value. The excess of the Hall term above this critical 
value is then included via a new scheme for pure Hall diffusion. 
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1 INTRODUCTION 

Dynamically important magnetic fields are commonplace in 
astrophysics. In many cases, where these fields interact with 
fluids, researchers have assumed that the equations of ideal 
magnetohydrodynamics (MHD) are sufficient in modeling 
the evolution of the magnetic fields and the fiuids with which 
they interact. There are clear examples, however, where the 
assumptions underpinning the equations of ideal MHD are 
not vahd. In dense molecular clouds, for example, the den- 
sity of charged p articles can be much low er than that of the 
neutral species JCiolek fc Robergell200l hereafter CR02). 
Under these conditions, coupling between the motions of 
the fiuids and the magnetic field is not perfect, and dif- 
fusive effects become significant. Similarly, ideal MHD is 
not believed to be valid in accretion disks around young 
stellar objects JWardldl200i) . The latter point is particu- 
larly interesting given the importance attached to the in- 
teraction between accretion disks and mag netic fields in the 
launching of stellar jets and outflows (e.g. IShu et all 
Fendt fc Camenzindl [l995 lOuved. Pudritz fc Stond 
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Ler^^tH ' ^Qgr ' ^^reiral l2004l) . When modeling systems 
such as these therefore, a full multifluid treatment permit- 
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tur- 



ting relative motions between different component species 
should be adopted. 

Many authors llTothI Il9 94": 'Smith & Mac Low' Il997l : 
IStonj|l997l : IChieze. Pineau des Forets fc FloweE.,.199A have 
suggested schemes for numerically integrating the multi- 
fluid equations in the limit of pure ambipolar diffusion. In 
this regime the charged species are firmly tied to the mag- 
netic field lines as they diffuse through the neutral gas. 
The problem becomes more technically challenging, how- 
ever, when charged species may be loosely attached to the 
field lines and Hall diffusion can become important. No- 
tably, it is thought that Hall diffusion may play an im- 
portant rol e in environments such as t he surfaces of neu- 
tron stars ifH ollerbach fc Riidigeil 120041) . protostellar disks 
iWard le 2004), and dense molecular clouds (CR02). 

In their numerical studies of molecular clouds, CR02 as- 
sumed that the ionisation fraction is low and that the inertia 
of the charged particles may be neglected. They were then 
able to integrate the governing equations for a multifiuid 
problem including the presence of several sp ecies of charge- 
carrymg gram. Separately, Sano fc Stone (l2002al : l2002bl) 
performed multifluid calculations designed to examine the 
Hall effect in the context of the magnetorotational instabil- 
ity in accretion disks. However, both of the schemes used 
by these authors are subject to a rather stringent stability 
criterion which requires that t he timestep tends to zero as 
the HaU effect becomes large llFalld l2003l hereafter F03). 
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To circumvent this constraint F03 presents a scheme em- 
ploying an implicit method of integrating the magnetic field 
equation. This has the advantage of allowing timesteps up 
to the limit dictated by the hyperbolic components of the 
equations. However, since large scale multifluid simulations 
are of obvious interest, the inherent difficulty of parallelising 
implicit schemes becomes a serious disadvantage. 

In this work we present a fully explicit numerical scheme 
for solving the multifiuid equations describing a weakly 
ionised plasma. The usual stability restrictions are relaxed 
through a combi nation of a technique known as Supe r Time 
Stepping (STS) iAlexiades. Amiez fc Gremaudlll996h and a 
new method which we call the Hall Diffusion Scheme (HDS). 
Crucially, since the scheme is explicit, it is straightforward 
to parallelise and to implement on top of an adaptive mesh 
refinement (AMR) engine. 

In Section|21the governing equations are described; Sec- 
tion |3| contains a detailed description and analysis of the nu- 
merical scheme; Section |1| contains numerical tests demon- 
strating the reliability of the scheme; and in Section |S] the 
relevance of this work is discussed. 



2 THE MULTIFLUID EQUATIONS 

We assume a weakly ionised plasma such that the mass 
density is dominated by the neutral component of the gas. 
Then, relative to the scale-length of the system, if particles 
of a given charged species have small mean free paths in the 
neutral gas, or small Larmor radii, their pressure and inertia 
may be neglected. 

For convenience it is assumed there is no mass trans- 
fer between species. It is straightforward, however, to insert 
the necessary terms for a more general treatment to include 
mass transfer (for example, see F03 and CR02) if desired. 
The equations governing the evolution of the multifluid sys- 
tem (CR02; F03) can then be written as 



9pi d 



J X B, 



(1) 



(2) 



dei d 
dt dx 



J-E + Y,H^, (3) 



dB dM _ d ^dB 

dt dx dx dx ' 



aipi {E + q- X B) + pipiKi i (Qi - = 0, 



Hi + Gii + OipiQi ■ E ^0, 



dB^ 
dx 



= 0, 



(4) 



(5) 



(6) 



(7) 



^ aiPi 



= 0, 



^ aiPiQ^ = J. 



(8) 



(9) 



The subscripts denote the species, with a subscript of 1 indi- 
cating the neutral fluid. The variables pi, q- = {ui,Vi,Wi)'^ 
and Pi are the mass density, velocity and pressure of 
species i. The identity matrix, current density and magnetic 
flux density are represented by /, J, B respectively. Ki i de- 
scribes the coUisional interaction between species i and the 
neutral fluid, ai is the charge-to-mass ratio for species i, Gn 
is the energy transfer rate from species i to the neutral fluid. 
Hi is the energy source or sink appropriate to species i, R 
is the resistivity matrix and M is the hyperbolic flux of B. 
See F03 and CR02 for a more detailed description of these 
terms. Note that in general Ki i and d i may depend on 
the temperatures and relative velocities of the interacting 
species. Equations Q to Q are the equations governing 
the conservation of mass, neutral momentum, neutral en- 
ergy, magnetic flux, charged momentum, and charged en- 
ergy. Equations (O to @ describe the divergence of B, 
charge neutrality, and current respectively. 

From Faraday's law in one dimension dBx/dt = so 
that the trivial Bx component may be dropped from equa- 
tion 1^. The hyperbolic flux is then 



M = {uiBy - viBx,U]_Bz - wiBx 
and the resistivity matrix is 



(10) 



R 



{ro - "Ta)^ +rA 



, .ByB, , x-By , ]^ ' 

[rA-ro)-^ rn^ (ro~rA)^ + rA J 



where ro, ru and ta are the Ohmic, Hall and ambipolar 
resistivities respectively and are defined by 



ro = 



TA 



1 



OA 



'H 



v2 ' 



with conductivities 

JV 



OA = 



OiiptPi, 

i = 2 

i_ \ - OipiPi 

B 



+ 0? 



where the Hall parameter for species i is given by 

a^B 



Kiipi 



(12) 
(13) 
(14) 

(15) 
(16) 
(17) 

(18) 
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3 NUMERICAL APPROACH 



3.2 Magnetic diffusion 



3.1 The gas equations 

Assuming a piecewise constant solution at time on a 
uniform mesli of spacing h, the solution at a later time 
f""*"^ = + r is sought. The state in cell j represents the 
volume average over (j — 1/2) /i < a; < (j + l/2)/i. 

It should first be noted that the charged particle pres- 
sures'^ and velocities {p'l^^ and q^'^^ for i > 1) can be ob- 
tained algebraically through equations Q and JSJ. This pro- 
cedure is described in Appendix 1X1 

To obtain the full solution at time i^^^ , finite volume 
methods are applied to equations Q to Q. The time inte- 
gration is multiplicatively operator split into five operations, 
with each carried out to second order accuracy in space and 
time. The order is permuted over successive timesteps such 
that seco nd order temp oral accuracy is maintained over the 
full step (IStranelll968^ . In the following the five necessary 
operations for finite volume integration are described. 

(i) Equations Q to (with i = 1 for the mass equation) 
form a system of equations for the neutral gas. Working in 
terms of the primitive variables P = (pi, <7i, pi)"^, fiuxes 
are evaluated from a piecewise constant solution P" via 
a hydrodynamic Riemann solver. A time centred solution, 

obtained from these fluxes is then reconstructed to 
a second order piecewise-linear solution, p"+i/^^ using Van 
Albada nonlinear averaging for the gradients. Fluxes may 
then be derived from p^+^Z^ which are second order accu- 
rate insmce and time (for further details see, for example, 
lFaIlelll99J) . These fluxes are then applied to the conserved 
variables. 

(ii) The source terms on the right hand sides of equa- 
tions ^ and ^ are applied. 

(iii) The charged particle mass fluxes are applied using 
equation ^ with i > 1 in a second order upwind procedure 
similar to that used for the neutral gas. 

(iv) The hyperbolic flux on the left hand side of equa- 
tion ^ is applied via a centred approximation 



M 



3+1/2 = 2 (-'^J+i 



(19) 



Splitting the hyperbolic flux term dM /dx from the induc- 
tion equation |0J and linearising yields 



dB 
'dt 



R 



(20) 



Note that the linearised form is assumed for convenience 
in the following analysis and in practise generalised discreti- 
sations of the nonlinear diffusion term are used. 



3.2.1 Standard discretisation 

The usual explicit discretisation for a diffusion term applied 
to equation II2UI I yields 



R';{B]+,-2B] +B]_^) 



(21) 



Assuming ro to be negligible, the relative impor- 
tance of the remaining resistivities can be parameterised by 
ri = rA/\rH\. F03 showed the above scheme has an amplifi- 
cation matrix with eigenvalues which are real when rj > r]* 
and complex otherwise. The transition point 77* is given by 



77* = 2lcos0|/sin^6l 



(22) 



where 6 is the pitch angle of the field with respect to the 
i-axis. In the real regime, the stability limit on the timestep 



2vTT^ 



r;(l -I- coa'^e) + 2|cos6'| v'(r?/j?*)^ - 1 



(23) 



where r = r/r"'" and is the characteristic cell crossing 
time for diffusion perpendicular to the magnetic field given 
by 



(24) 



2\rH\^/TTv^ 

However, below the transition point the stability limit be- 



This has the disadvantage of not coupling the bulk fiuid 
to the magnetic field through a Riemann problem, however, 
it is necessary in order that purely hydrodynamic subshocks 
may be properly captured. As remarked by F03, as long 
as the magnetic field appears continuous on the grid, as 
should be the case with finite resistivities, this is perfectly 
acceptable. 

(v) The resistive term on the right hand side of equa- 
tion @ is applied. Discussion of this procedure is deferred 
to the following section since it is of special interest. 



1 + cos'' 



(25) 



In either case the stable timestep limit goes as h?' since this 
is an explicit discretisation of a diffusion equation, however, 
a potentially more severe constraint is that while this limit 
increases as r) ~* rf in the real regime, it rapidly drops to 
zero as r; — > in the complex regime. 

3.2.2 Numerical Strategy 

Our strategy is to split rn into two parts such that 



^ It is actually the charged species' temperatures which are de- 
rived as their pressures are not explicitly necessary under the 
assumptions made here. 



rn = rH + 



(26) 



where = ;p»'£f is the maximum allowable Hall resistiv- 
ity in the real regime and is the excess. The induction 
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equation is then integrated in two parts using a technique 
to accelerate the timestepping for the standard discretisa- 
tion with Hall resistivity r"^. The excess Hall resistivity r\j 
is then applied using a different discretisation with suitable 
stability properties. 



3.2.3 Super Time Stepping 

STS is a technique which can be used to accelerate explicit 
schemes for parabolic problems. Essentially a Runge-Kutta- 
Chebyshov method, it has been known for some time (see 
L^-lcxiadcs, Amicz & Grcmaud 1996), although it remains 
relatively unknown in computational astrophysics. 

A superstep r^^^ is a composite timestep built up from 
a series of A^sts substeps such that 



where the cosine term is absorbed by defining d\i = r\jcos0. 
It seems to make little difference which component is ad- 
vanced first. 

For clarity of notation the superscript b is dropped 
from the following analysis of the stability properties of the 
scheme. The resistance matrix for pure Hall diffusion is 



R 











Assuming a numerical wave of the form 



(31) 



(32) 



in equations 11291 and Il3m yields an amplification matrix 



A = 



1 

dn 



-4 



(33) 



iVgrps 



where du = (.dn and 



(27) 



Judicious choice of the dr, yields stability for the super- 
step while the normal stability restrictions on the individual 
substeps are relaxed. Exploiting the properties of Chebyshev 
polynomials provides a set of optimal values for the substeps 
given by 



dri 



-1 + u) cos 



2 j - 1 TT 

AsTS 2 



+ 1 + 1^ 



(28) 



where is the normal explicit timestep limit and u is a 
damping factor. Note that r^"^^ ^ A'Itst^ as ly ^ 0. The 
method is uns table in the limit v = 0. For a more detailed 
discussion, see lAlexiades. Amiez fc GremaudI l|l99(fl and ref- 
erences therein. 

In order to apply STS to second order in time Richard- 
son extrapolation in used. 



3.2.4 Hall Diffusion Scheme 

Having advanced the induction equation with a Hall resistiv- 
ity rfi, it is necessary to find an efficient scheme to impose 
the excess Hall diffusion r'^j. Since multiplicative operator 
splitting yields a composite scheme with an amplification 
factor equal to the product of the amplification factors of 
the basis schemes this task can be reduced to one of finding 
a scheme for pure Hall diffusion. 

The key observation to make is that R has zero entries 
on the diagonal when pure Hall diffusion is being considered. 
With this in mind, equation 12111 may be used to advance 
one component of the magnetic field explicitly, followed by 
an implicit-like discretisation of the alternate component. 
We call this the Hall Diffusion Scheme (HDS) as we are 
not aware of an instance of this approach elsewhere in the 
literature. Hence the discretisation of equation 12111 for the 
pure Hall excess becomes 



2b: 



followed by 



,n + l 



/!,2 



b: 



-^yj-i) 



(29) 



(30) 



2r(l — cosoj) 



/l2 



The eigenvalues of A are given by 

.1 

'2 



A = 1 - idlf ± i^dH\/'i - d'i 



(34) 



(35) 



and hence HDS is neutrally stable for |iijf | < 2. Taking the 
most restrictive case of u — n gives a stable timestep limit 
of 



\cose\{l-rj/v*)' 



(36) 



Note that f™^ l/|cosS| as 77 in contrast to the 
standard discretisation for which f — > 0. 

The extension of the HDS to more than one dimension is 
straightforward although we defer a detailed discussion to a 
later paper. For an outline of the scheme in three dimensions 
the reader is referred to Appendix IHl 

In practice, ordinary (unaccelerated) subcycling of 
HDS, using Ahds subcycles, is applied in conjunction with 
STS. This compound scheme (referred to as "STS/HDS" 
hereafter) usually allows the timestep limit imposed by 
the hyperbolic terms to be reached efficiently (see Subsec- 
tion 



4 NUMERICAL TESTS 

Following F03, the dynamic algorithm described here is 
tested against solutions of the steady isothermal multifluid 
equations. These steady state equations are solved using an 
independent code, the details of which are outlined in Ap- 
pendix |n| The conditions for each of the tests are given in 
Table El 



4.1 Case A: Ambipolar Dominated 



In this test ro = 2 x 10"'^, rn = 1.16 x 10"'' and ta = 0.068 
giving r] = 5.86 x 10"^ and hence it can be expected that 
ambipolar diffusion will dominate the solution. Fig.^shows 
plots of the X component of the neutral velocity, along with 
By for both the dynamic and steady state solutions. The 
calculation shown has h = 5 x 10~^ . It can be seen that the 
agreement between the two solutions is extremely good. 
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Case A 


















Right State 


Pi 


= 1 


<Ji = (—1.751, 0, 0) 


B = (1, 0.6, 0) 


P2 = 


5 X 10 


P3 


= 1 X 10 


Left State 


Pi 


= 1.7942 


Qi = (—0.9759,-0.0561,0) 


B = (1, 1.74885, 0) 


P2 = 


8.9712 X 10 


P3 


= 1.7942 X 10 




02 


= —2 X lO"""^ 


as = 1 X 10* 


K21 = 4 X 10^ 




= 2 X 10* 


a - 


= 0.1 




V - 


= 0.05 


JVSTS = 5 


A^HDS = 










Case B 


















Right State 


As 


case A 














Lett btate 


As 


case A 
















OL2 


= — 2 X 10^ 


as = 1 X 10° 


K21 = 4 X lO'' 


K31 


= 2.5 X 10° 


a - 


= 0.1 




V - 


= 


^STS = 1 












Case C 


















Right State 


Pi 


= 1 


= (-6.7202,0,0) 


B = (1,0.6,0) 


P2 = 


5 X 10-8 


P3 


= 1 X 10-3 


Left State 


Pi 


= 10.421 


= (-0.6449, -1.0934,0) 


B= (1,7.9481,0) 


P2 = 


5.2104 X 10-^ 


p-i 


= 1.0421 X 10-2 




02 


= -2 X 10^2 


as = 1 X 10* 


iCa 1 = 4 X 105 


X31 


= 2 X 10* 


a - 


= 1 




V - 


= 0.05 


AfSTS = 15 


A^HDS = 











Table 1. Test calculation parameters. 



Since the algorithm is designed to be second order it is 
worthwhile measuring the convergence rate of the dynamic 
solution against the solution from the steady state solver. 
The comparison is made by minimising the LI error norm, 
ei, between a section of the dynamical solution and the 
steady state solution. Working from the downstream side, 
the section xl < x < xr is fixed about the point x* where 
the deviation from the downstream state first exceeds 50% of 
the maximum variation in the solution. Using xl = a;* — 0.44 
and XR = a;* +0.56 yields ei = 3.90 x lO"'^ for /i = 5 x 10"^ 
and ei = 1.56 x 10-* for /i = 1 x 10-^. This gives ei oc h? ° 
- showing second order convergence as expected. 



4.2 Case B: Hall Dominated 

The Hall term dominates in this test, requiring the Hall 
diffusion to be split and applied in part via HDS. The pa- 
rameters are ro = 2 x 10"®, rn = 0.0116, ta = 5.44 x 10"* 
with r) — 0.0046 ^ 1? Fig.|5|shows the results of the calcu- 
lations for the test with h — 2 x 10"^. For standard explicit 
codes the conditions lead to prohibitive restrictions on the 
timestep. However, the use of HDS allows us to maintain a 
timestep close to the Courant limit imposed by the hyper- 
bolic terms throughout the calculations. 

As with case A, the dynamic solution is tested to ensure 
it has the correct second order convergence characteristics. 
With XL ~ X* — 0.15 and xr = x* + 0.95, we find ei = 
4.95 X 10-3 for /i = 2 X lO^^ and ei = 1.15 x lO^^ for 
h = 1 X 10""^ , giving ei oc h^'^ . Again, this is close to the 
second order convergence rate expected. 



4.3 Case C: Neutral subshock 

This test is similar to case A, but with a higher soundspeed 
and upstream fast Mach number. As a result, a subshock 
develops in the neutral fiow because the interactions between 
the charged particles and the neutrals are not strong enough 
to completely smooth out the strong initial discontinuity in 

^ If the Hall diffusion is increased much further, it appears that 
the approximation of negligible charged particle inertia breaks 
down. 




0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 

X 

Figure 1. Neutral fluid x- velocity and ^/-component of magnetic 
field for case A with h = 5 X 10~^. The solution from the steady 
state equations, as a line, is overplotted with points from the 
dynamic code. 



the neutral flow. The ability of the algorithm described to 
deal with discontinuities in the solution is therefore tested. 

Fig. |3| shows the results of the calculations for h = 
1 X 10"'^. The subshock in the neutral flow is clearly visible 
as a discontinuity in ui , while there is no corresponding dis- 
continuity in By. Fig. |1| contains a plot of the x component 
of the velocity of the negatively charged fluid. As expected, 
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u +1 t+ t| t 

t V' ^ * ' 



Dynamic solution 
Steady solution 



Ul 




Dynamic solution 
Steady solution 





Dynamic solution 
Steady solution 



Figure 2. Neutral fluid a;-velocity and j/-component of magnetic 
field for case B with h = 2 x 10~^. The solution from the steady 
state equations, as a line, is overplotted with points from the 
dynamic code. 



Figure 3. Neutral fluid x-velocity and ^/-component of magnetic 
field for case C with h = 1 X 10~^. The solution from the steady 
state equations, as a line, is overplotted with points from the 
dynamic code. 



there is no discontinuity in this variable, but there are some 
oscillations at the point where the discontinuity in the neu- 
tral flow occurs. These errors are remarkably similar to those 
encountered by F03 and do not alTect the global solution. 

It can be expected that, since there is a discontinuity in 
the solution of this test and a MUSCL-type approach is used, 
the rate of convergence of the dynamic solution will be close 
to first order, at least for resolutions high enough to discern 
the subshock in the solution. In this test xl = x* —0.13 and 
XR = x* + 0.15. We find ei = 3.41 x 10"^ for /i = 5 x 10~^ 
and ei = 5.25 x 10"^ for /i = 1 x 10"^ yielding ei oc h^'^^ - 
close to the first order expected, although clearly the error 
from around the subshock is not completely dominating at 
this resolution. At h = 5 x 10"" we find ei = 2.73 x 10"^ 
giving ei oc h"'^* with respect to the error at h — 1 x 10"'^. 
We suspect that the deviation from first order is due to a 
discontinuity in the electric field at the subshock causing 
an error in the charged velocities since smoothing the solu- 
tion with some artificial viscosity is found to improve the 
convergence. 



U2 




Dynamic solution 
Steady solution 



X 

Figure 4. Negatively charged fiuid x-velocity for case C with 
?i = 1 X 10"'^. The solution from the steady state equations, as a 
line, is overplotted with points from the dynamic code. 
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Case A 


Case B 


Case C 


STS/HDS 


1.9 


14.8 


1.9 


Implicit 


1.9 


23.3 


2.7 



Table 2. The speed-up factors in CPU time usage achieved 
via the impUcit and STS/HDS discretisations of the induction 
equation relative to the standard explicit discretisation. 

4.4 Comparative timings 

In this section comparison is made between the perfor- 
mances of standard explicit, STS/HDS, and implicit (Crank- 
Nicolson) discretisations of the induction equation. The dif- 
ferent methods are applied in otherwise identical codes to 
the iiigli resolution trials of the preceding test cases. Since 
the neutral gas equations are treated explicitly in all cases, 
the corresponding Courant condition on the integration of 
the hyperbolic terms imposes a hard limit on the timestep. 

As a benchmark, we use the standard explicit discretisa- 
tion subcycled to the same degree as the STS /HDS method. 
The speed-up factors of the STS/HDS and implicit meth- 
ods in terms of CPU time usage are presented in Table |21 
Clearly, either technique offers a significant improvement in 
efficiency and both achieve timesteps close to the limit in- 
troduced by the hyperbolic terms. The implicit method is 
slightly faster for case C due to the iiigii degree of subcy- 
cling used for the STS and significantly so for case B because 
of the very large Hall term. Otherwise, the STS/HDS and 
implicit methods yield similar speed-up factors indicating 
that overall efficiency is dominated by the other parts of 
the schemes. It should be emphasised that these are steady 
state problems which suit implicit methods particularly well 
and for non-steady state problems accuracy constraints may 
reduce the efficiency of implicit schemes. 



5 CONCLUSIONS 

A new explicit scheme for integrating the multifluid equa- 
tions in the limit of low ionisation has been presented. The 
usual explicit stability limit imposed by the induction equa- 
tion is relaxed by means of the Super Time Stepping al- 
gorithm applied for a portion of the Hall diffusion up to a 
critical limiting value. 

Beyond this limiting value the standard explicit dis- 
cretisation becomes subject to a stability constraint requir- 
ing that the timestep vanish as the Hall diffusion becomes 
large. In order to circumvent this constraint, the excess Hall 
diffusion above the critical value is split off and applied via a 
new method which we have called the Hall Diffusion Scheme. 

It has been demonstrated that, for the case of an 
isothermal flow, the algorithm is accurate and converges to 
second order when the solution is smooth and to flrst order 
when the solution contains a discontinuity. The extension 
of this scheme to non-isothermal flow does not present any 
obvious difficulties, although a modiflcation of the discreti- 
sation used for the magnetic flux evolution is necessary. 

Since all discretisations used in the scheme presented 
here are explicit, it is a straightforward matter to implement 
in a multidimensional parallelised codes using adaptive mesh 
refinement. This is a crucial advantage for large-scale sim- 



ulations of astrophysical systems in which multifluid effects 
are thought to be important such as dense molecular clouds 
and protostellar accretion disks. 
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APPENDIX A: CHARGED VELOCITIES 

For this work the coUisional coefficients Ki i are assumed 
to be independent of velocities and temperatures. The fol- 
lowing derivation (S.A.E.G. Falle, private communication) 
is included for completeness. 

Transforming to the frame co-moving with the neutral 
gas, equation @ can be written as 

q[xB~ K^q[ = E' q^xB (Al) 

where Ki = p\Ki i/ai and E' = E + q^ x B. 

Then choosing i = 2 as a reference species, the general 
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solutions for the remaining charged species' velocities are 
given by 



where 




-By 

Bx 



(A2) 



(A3) 



To derive the charged velocities, all that remains is for 
the reference velocity to be evaluated, this can be done by 
using equation Q and Ampere's law to give 



By 

B"+' 



— Bx ~\'T{GxyBy 

= By + T{GyzB" 

= B" + t{G"xBx 



+ G" 



zB"), 

z^x ) ^ 



(B5) 
(B6) 
(B7) 



where is the discretised form of the operator G at time 
level n. 

The generalised HDS scheme in three dimensions is 
analogous in construction to the one-dimensional case in 
that equation l|B5fl is an explicit first step and equation ljB7^ 
is an implicit-like final step. Additionally, we now have an 
intermediate step of mixed explicit /implicit character. Nu- 
merical tests indicate that the method retains its favourable 
stability properties in three dimensions. 



Q2 



E 



Clip 



a--2.pi 



-a: 



V xB 

a2p2 



(A4) 



If the collisional coefficients are in fact dependent on 
the velocities of the charged species, this procedure can be 
carried out iteratively using the values from the previous 
timestep as a starting point. 

Should the collisional coefficients depend on the temper- 
atures of the charged species, some additional calculation is 
necessary before the next iteration: using equation © and 
inserting the specific form of the function Gn, N ~ 1 equa- 
tions are obtained which may be solved readily for the A'^ — 1 
charged temperatures. 

Finally, it is worth noting that superior results are ob- 
tained by interpolating the primitive quantities to the cell 
edges before calculating the charged velocities rather than 
by calculating the velocities at the cell centres and interpo- 
lating from these to the edges. 



APPENDIX B: HDS IN THREE DIMENSIONS 

Equation JKJ can be used in conjunction with with equa- 
tion ^ to write the electric field for pure Hall diffusion as 



E^rn 



J X B 

B 



Then, using Faraday's law, we can write 
dB 



dt 



-rj/V X {J X b) 



(Bl) 



(B2) 



where b = B/B. 

This equation can be expanded out and linearised to give 



dB 

dt 



GB 



(B3) 



where, using J = V x B, the matrix G is given by 



APPENDIX C: STEADY-STATE SOLVER 

Assuming an isothermal flow, as is the case for the tests 
presented in this work, setting all derivatives with respect 
to time to zero in the multifluid equations gives us 



pi Ml + pi + 



piUi 

B^ 



— Qi = constant, 
= Px = constant, 



piUivi — BxBy — Py = constant, 



PlUlWl 



BxB, — Pz ~ constant. 



(CI) 

(C2) 

(C3) 
(C4) 



In addition the reduced momentum equations for the 
charged species (equation Q) yield 3 equations for each 
charged species, and the charge neutrality condition is also 
used. Finally, the equation for B yields 

.dB 



M - Mr = R 



dx 



(C5) 



with M-R (= M-l) being the flux in the right (left) state. 

For the cases considered here, with two charged species, 
the above equations constitute one ordinary differential 
equation for B and seven equations which, once B is known 
at a given point in space, can be used to solve for all the 
other variables. The ODE for B is solved using the Runge- 
Kutta method of order 4. 

The initial conditions (at a; = 0) are a saddle point of 
the ODE for B. These conditions are perturbed slightly and 
the system then evolves through phase space to a sink point. 



G = -rH{b-V)V X 



(B4) 



Hence G is antisymmetric and we can write the generalised 
HDS scheme as 
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